1 /*
2 * Copyright (C) 2011 The Guava Authors
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16
17 package com.google.common.math;
18
19 import static com.google.common.base.Preconditions.checkArgument;
20 import static com.google.common.math.MathPreconditions.checkNonNegative;
21 import static com.google.common.math.MathPreconditions.checkPositive;
22 import static com.google.common.math.MathPreconditions.checkRoundingUnnecessary;
23 import static java.lang.Math.min;
24 import static java.math.RoundingMode.HALF_EVEN;
25 import static java.math.RoundingMode.HALF_UP;
26
27 import com.google.common.annotations.GwtCompatible;
28 import com.google.common.annotations.VisibleForTesting;
29
30 import java.math.RoundingMode;
31
32 /**
33 * A class for arithmetic on values of type {@code long}. Where possible, methods are defined and
34 * named analogously to their {@code BigInteger} counterparts.
35 *
36 * <p>The implementations of many methods in this class are based on material from Henry S. Warren,
37 * Jr.'s <i>Hacker's Delight</i>, (Addison Wesley, 2002).
38 *
39 * <p>Similar functionality for {@code int} and for {@link BigInteger} can be found in
40 * {@link IntMath} and {@link BigIntegerMath} respectively. For other common operations on
41 * {@code long} values, see {@link com.google.common.primitives.Longs}.
42 *
43 * @author Louis Wasserman
44 * @since 11.0
45 */
46 @GwtCompatible(emulated = true)
47 public final class LongMath {
48 // NOTE: Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
49
50 /**
51 * Returns {@code true} if {@code x} represents a power of two.
52 *
53 * <p>This differs from {@code Long.bitCount(x) == 1}, because
54 * {@code Long.bitCount(Long.MIN_VALUE) == 1}, but {@link Long#MIN_VALUE} is not a power of two.
55 */
56 public static boolean isPowerOfTwo(long x) {
57 return x > 0 & (x & (x - 1)) == 0;
58 }
59
60 /**
61 * Returns 1 if {@code x < y} as unsigned longs, and 0 otherwise. Assumes that x - y fits into a
62 * signed long. The implementation is branch-free, and benchmarks suggest it is measurably
63 * faster than the straightforward ternary expression.
64 */
65 @VisibleForTesting
66 static int lessThanBranchFree(long x, long y) {
67 // Returns the sign bit of x - y.
68 return (int) (~~(x - y) >>> (Long.SIZE - 1));
69 }
70
71 /**
72 * Returns the base-2 logarithm of {@code x}, rounded according to the specified rounding mode.
73 *
74 * @throws IllegalArgumentException if {@code x <= 0}
75 * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x}
76 * is not a power of two
77 */
78 @SuppressWarnings("fallthrough")
79 // TODO(kevinb): remove after this warning is disabled globally
80 public static int log2(long x, RoundingMode mode) {
81 checkPositive("x", x);
82 switch (mode) {
83 case UNNECESSARY:
84 checkRoundingUnnecessary(isPowerOfTwo(x));
85 // fall through
86 case DOWN:
87 case FLOOR:
88 return (Long.SIZE - 1) - Long.numberOfLeadingZeros(x);
89
90 case UP:
91 case CEILING:
92 return Long.SIZE - Long.numberOfLeadingZeros(x - 1);
93
94 case HALF_DOWN:
95 case HALF_UP:
96 case HALF_EVEN:
97 // Since sqrt(2) is irrational, log2(x) - logFloor cannot be exactly 0.5
98 int leadingZeros = Long.numberOfLeadingZeros(x);
99 long cmp = MAX_POWER_OF_SQRT2_UNSIGNED >>> leadingZeros;
100 // floor(2^(logFloor + 0.5))
101 int logFloor = (Long.SIZE - 1) - leadingZeros;
102 return logFloor + lessThanBranchFree(cmp, x);
103
104 default:
105 throw new AssertionError("impossible");
106 }
107 }
108
109 /** The biggest half power of two that fits into an unsigned long */
110 @VisibleForTesting static final long MAX_POWER_OF_SQRT2_UNSIGNED = 0xB504F333F9DE6484L;
111
112 // maxLog10ForLeadingZeros[i] == floor(log10(2^(Long.SIZE - i)))
113 @VisibleForTesting static final byte[] maxLog10ForLeadingZeros = {
114 19, 18, 18, 18, 18, 17, 17, 17, 16, 16, 16, 15, 15, 15, 15, 14, 14, 14, 13, 13, 13, 12, 12,
115 12, 12, 11, 11, 11, 10, 10, 10, 9, 9, 9, 9, 8, 8, 8, 7, 7, 7, 6, 6, 6, 6, 5, 5, 5, 4, 4, 4,
116 3, 3, 3, 3, 2, 2, 2, 1, 1, 1, 0, 0, 0 };
117
118 // halfPowersOf10[i] = largest long less than 10^(i + 0.5)
119
120 /**
121 * Returns the greatest common divisor of {@code a, b}. Returns {@code 0} if
122 * {@code a == 0 && b == 0}.
123 *
124 * @throws IllegalArgumentException if {@code a < 0} or {@code b < 0}
125 */
126 public static long gcd(long a, long b) {
127 /*
128 * The reason we require both arguments to be >= 0 is because otherwise, what do you return on
129 * gcd(0, Long.MIN_VALUE)? BigInteger.gcd would return positive 2^63, but positive 2^63 isn't
130 * an int.
131 */
132 checkNonNegative("a", a);
133 checkNonNegative("b", b);
134 if (a == 0) {
135 // 0 % b == 0, so b divides a, but the converse doesn't hold.
136 // BigInteger.gcd is consistent with this decision.
137 return b;
138 } else if (b == 0) {
139 return a; // similar logic
140 }
141 /*
142 * Uses the binary GCD algorithm; see http://en.wikipedia.org/wiki/Binary_GCD_algorithm.
143 * This is >60% faster than the Euclidean algorithm in benchmarks.
144 */
145 int aTwos = Long.numberOfTrailingZeros(a);
146 a >>= aTwos; // divide out all 2s
147 int bTwos = Long.numberOfTrailingZeros(b);
148 b >>= bTwos; // divide out all 2s
149 while (a != b) { // both a, b are odd
150 // The key to the binary GCD algorithm is as follows:
151 // Both a and b are odd. Assume a > b; then gcd(a - b, b) = gcd(a, b).
152 // But in gcd(a - b, b), a - b is even and b is odd, so we can divide out powers of two.
153
154 // We bend over backwards to avoid branching, adapting a technique from
155 // http://graphics.stanford.edu/~seander/bithacks.html#IntegerMinOrMax
156
157 long delta = a - b; // can't overflow, since a and b are nonnegative
158
159 long minDeltaOrZero = delta & (delta >> (Long.SIZE - 1));
160 // equivalent to Math.min(delta, 0)
161
162 a = delta - minDeltaOrZero - minDeltaOrZero; // sets a to Math.abs(a - b)
163 // a is now nonnegative and even
164
165 b += minDeltaOrZero; // sets b to min(old a, b)
166 a >>= Long.numberOfTrailingZeros(a); // divide out all 2s, since 2 doesn't divide b
167 }
168 return a << min(aTwos, bTwos);
169 }
170
171 @VisibleForTesting static final long FLOOR_SQRT_MAX_LONG = 3037000499L;
172
173 static final long[] factorials = {
174 1L,
175 1L,
176 1L * 2,
177 1L * 2 * 3,
178 1L * 2 * 3 * 4,
179 1L * 2 * 3 * 4 * 5,
180 1L * 2 * 3 * 4 * 5 * 6,
181 1L * 2 * 3 * 4 * 5 * 6 * 7,
182 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8,
183 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9,
184 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10,
185 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11,
186 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12,
187 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13,
188 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14,
189 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15,
190 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16,
191 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17,
192 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18,
193 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18 * 19,
194 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18 * 19 * 20
195 };
196
197 /**
198 * Returns {@code n} choose {@code k}, also known as the binomial coefficient of {@code n} and
199 * {@code k}, or {@link Long#MAX_VALUE} if the result does not fit in a {@code long}.
200 *
201 * @throws IllegalArgumentException if {@code n < 0}, {@code k < 0}, or {@code k > n}
202 */
203 public static long binomial(int n, int k) {
204 checkNonNegative("n", n);
205 checkNonNegative("k", k);
206 checkArgument(k <= n, "k (%s) > n (%s)", k, n);
207 if (k > (n >> 1)) {
208 k = n - k;
209 }
210 switch (k) {
211 case 0:
212 return 1;
213 case 1:
214 return n;
215 default:
216 if (n < factorials.length) {
217 return factorials[n] / (factorials[k] * factorials[n - k]);
218 } else if (k >= biggestBinomials.length || n > biggestBinomials[k]) {
219 return Long.MAX_VALUE;
220 } else if (k < biggestSimpleBinomials.length && n <= biggestSimpleBinomials[k]) {
221 // guaranteed not to overflow
222 long result = n--;
223 for (int i = 2; i <= k; n--, i++) {
224 result *= n;
225 result /= i;
226 }
227 return result;
228 } else {
229 int nBits = LongMath.log2(n, RoundingMode.CEILING);
230
231 long result = 1;
232 long numerator = n--;
233 long denominator = 1;
234
235 int numeratorBits = nBits;
236 // This is an upper bound on log2(numerator, ceiling).
237
238 /*
239 * We want to do this in long math for speed, but want to avoid overflow. We adapt the
240 * technique previously used by BigIntegerMath: maintain separate numerator and
241 * denominator accumulators, multiplying the fraction into result when near overflow.
242 */
243 for (int i = 2; i <= k; i++, n--) {
244 if (numeratorBits + nBits < Long.SIZE - 1) {
245 // It's definitely safe to multiply into numerator and denominator.
246 numerator *= n;
247 denominator *= i;
248 numeratorBits += nBits;
249 } else {
250 // It might not be safe to multiply into numerator and denominator,
251 // so multiply (numerator / denominator) into result.
252 result = multiplyFraction(result, numerator, denominator);
253 numerator = n;
254 denominator = i;
255 numeratorBits = nBits;
256 }
257 }
258 return multiplyFraction(result, numerator, denominator);
259 }
260 }
261 }
262
263 /**
264 * Returns (x * numerator / denominator), which is assumed to come out to an integral value.
265 */
266 static long multiplyFraction(long x, long numerator, long denominator) {
267 if (x == 1) {
268 return numerator / denominator;
269 }
270 long commonDivisor = gcd(x, denominator);
271 x /= commonDivisor;
272 denominator /= commonDivisor;
273 // We know gcd(x, denominator) = 1, and x * numerator / denominator is exact,
274 // so denominator must be a divisor of numerator.
275 return x * (numerator / denominator);
276 }
277
278 /*
279 * binomial(biggestBinomials[k], k) fits in a long, but not
280 * binomial(biggestBinomials[k] + 1, k).
281 */
282 static final int[] biggestBinomials =
283 {Integer.MAX_VALUE, Integer.MAX_VALUE, Integer.MAX_VALUE, 3810779, 121977, 16175, 4337, 1733,
284 887, 534, 361, 265, 206, 169, 143, 125, 111, 101, 94, 88, 83, 79, 76, 74, 72, 70, 69, 68,
285 67, 67, 66, 66, 66, 66};
286
287 /*
288 * binomial(biggestSimpleBinomials[k], k) doesn't need to use the slower GCD-based impl,
289 * but binomial(biggestSimpleBinomials[k] + 1, k) does.
290 */
291 @VisibleForTesting static final int[] biggestSimpleBinomials =
292 {Integer.MAX_VALUE, Integer.MAX_VALUE, Integer.MAX_VALUE, 2642246, 86251, 11724, 3218, 1313,
293 684, 419, 287, 214, 169, 139, 119, 105, 95, 87, 81, 76, 73, 70, 68, 66, 64, 63, 62, 62,
294 61, 61, 61};
295 // These values were generated by using checkedMultiply to see when the simple multiply/divide
296 // algorithm would lead to an overflow.
297
298 static boolean fitsInInt(long x) {
299 return (int) x == x;
300 }
301
302 /**
303 * Returns the arithmetic mean of {@code x} and {@code y}, rounded toward
304 * negative infinity. This method is resilient to overflow.
305 *
306 * @since 14.0
307 */
308 public static long mean(long x, long y) {
309 // Efficient method for computing the arithmetic mean.
310 // The alternative (x + y) / 2 fails for large values.
311 // The alternative (x + y) >>> 1 fails for negative values.
312 return (x & y) + ((x ^ y) >> 1);
313 }
314
315 private LongMath() {}
316 }